Waves on an interface between two phase-separated Bose-Einstein condensates 
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called, correspondingly, weakly and strongly segregated 
phase. The former one takes place if 512 exceeds ,/<?i 1.922 
only slightly. In such a case the component interpene- 
tration depth is quite large and proportional to 7 ' , 
where 



7 = 312(511522) 1/2 - 1. 



I.E.Mazets 
Ioffe Physico-Technical Institute, 194021 St. Petersburg, Russia 

We examine waves localized near a boundary between two weakly segregated Bose-Einstein con- 
densates. In the case of a wavelength of order of or larger than the thickness of the overlap region 
the variational method is used. The dispersion laws for the two oscillation branches are found in 
analytic form. The opposite case of a wavelength much shorter than the healing length in the bulk 
condensate is also discussed. 



PACS numbers: 03.75.Fi, 05.30.Jp, 67.60.-g 

The studies of mixtures of two supcrfluid Bosc liquids 
ascend to the early work by Khalatnikov M who deter- 
mined the possible sound modes in such a system within a 
macroscopic (hydrodynamical) approach. Later a micro- 
scopic theory was developed: Bassichis pf considered a 
neutral mixture of two charged Bose gases with Coulomb 
interactions, and Nepomnyashchii |H] considered a system 
composed of bosons of two kinds interacting via short- 
range potentials. In Rcf. J3|, most closely related to the 
case of a two-component degenerate dilute atomic vapor, 
the two-branch spectrum of elementary excitations and 
following from it criterion for stability of the ground state 
were found. 

The advances in experiments on Bose-Einstein con- 
densation of trapped alkali atoms gave rise to an in- 
terest to this subject. The ground state configuration 
of a two-component atomic Bose-Einstein condensate 
(BEC) in the presence of a harmonic confining potential 
at zero temperature was calculated by Ho and Shenoy 
H in the Thomas-Fermi limit. The collective oscilla- 
tion frequencies for trapped binary BECs were also deter- 
mined |g-0]. Interesting numerical results are obtained 
by Pu and Bigelow M . These studies reveal that if the 
number of trapped atoms of each kind is large enough, 
the ground state physics can be qualitatively understood 
from the arguments valid for a case of absence of a 
trap fl. In the latter case, the ground state proper- 
ties are determined by a certain relation between the 
coupling constants. Namely, if the intercomponent re- 
pulsion (we do not consider attractive potentials in the 
present paper) is small enough, i.e., if gi2 < JgxiQii-, 
where 5^ = 2irh[mimj / (toj + irij)]~ 1 aij, aij is the corre- 
sponding s-wave scattering length of a pair of ultracold 
atoms when one of them is of the jth kind and another 
is of the ith kind, rrij is the mass of an atom for the jth 
component of the mixed BEC, ij — 1, 2, then the two 
degenerate Bose gases are miscible, i.e. they co-exist in 
all the volume. In the opposite case, 512 > -^511522, the 
two components are immiscible and separated in space. 
In the latter case, a new physics related to the intercom- 
ponent boundary arises. The steady-state energetics of 
such an interface was estimated by Timmermans H . Ao 
and Chui |10| performed more detailed analysis, in par- 
ticular, they found that there are two different regimes 
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The latter limiting case can be achieved provided that 
7 ^> 1, under such a condition the density profile at the 
intercomponent boundary has quite a different shape and 
its thickness is determined by individual healing lengths 
for the bulk condensates. 

There are also experimental works on two-component 
BECs made by the JILA and NIST group Jll|,|| and by 
the MIT group §§. The JILA and NIST group created 
binary mixtures of ultracold bosons by populating dif- 
ferent magnetic sublevels, \F, m,p) of 87 Rb. It is worth 
to note that a mixture of the |2. ,2) and |1, —1) exhibits 
behaviour of a miscible system UW (cf. the related the- 
oretical paper |L4J), and the mixture of the |2, 1) and 
|1, — 1) looks like an example of a weakly segregated sys- 
tem @. 

In the present paper, we consider normal modes for 
intercomponent boundary oscillations in the weakly seg- 
regated phase regime, i.e., for < 7 <C 1, since a study 
of excited states of an immiscible system was lacking up 
to now. The Bogoliubov spectrum of a translationally 
invariant mixed BEC was obtained in Ref. H, but it re- 
vealed only the fact that such a homogeneous system is 
dynamically unstable against long-wavelength perturba- 
tions and, hence, tends to become phase-separated. One 
can think that dynamics of a phase separated system 
with an interface between the components 1 and 2 in 
the long-wavelength regime is adequately described by 
a concept of a surface tension with the surface tension 
constant inferred from the ground state properties fl9|. 
However, a rigorous proof of this result was still lacking. 
In the present paper, we fill this gap and examine also 
the opposite (short-wavelength) regime. 

Let us introduce two macroscopic wave functions for 
the components j — 1,2 as y/fij<f>j(r,t)exp(—i[j,jt/h), 
where fij = hgjjUj is the chemical potential of the corre- 
sponding component, and rij is its bulk number density. 
In the present paper, we consider a case of absence of 



any external potential, so the coupled set of the Gross- 
Pitaevskii equations IfM is reduced to 



Since, by assumption, 7 is a small parameter, we use the 
following decomposition in series: 
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We consider a case of two immiscible BECs. We denote 
the ground state solution for <£j by 4>j(r)- Let the first 
component occupy a semi-infinite space left to the z = 
plane that can be associated with the intercomponent 
boundary, i.e., 4>\ — > 1 if z — > —00 and 0i — » if z — > 
+00. Respectively, 02 — ► if 2 — > — 00 and 02 — > 1 if 
z — * +00. Hydrostatic arguments lead to the equality of 
the bulk component pressures Q: 



1 2 1 2 
25n«i = 2322^2- 



(4) 



To get a qualitatively correct physical picture, it is 
sufficient to consider a symmetric case: mi = mi and 
.9n — 322- Then we can introduce the time unit (j}\\n\)~ x 
and the length unit yjh/{2migi\n\) and write Eqs.(|[ ||) 
in the dimensionless form (hereafter, for the sake of 
brevity, we do not introduce specific notation for the di- 
mensionless time and coordinate): 



i$i = -V 2 $i - $1 + |$iT *i + (1 + 7) |*a| $ i> ( 5 ) 

i$ 2 = -V 2 $ 2 ^ $2 



|$ 2 | 2 $2 + (l+7)l$i| 2 *2. (6) 



Obviously, the steady state solution of Eqs.(g, |6|) pos- 
sesses the following symmetry: 

<h(») = <h(-*)- (?) 

Now we have to determine the steady state order pa- 
rameters for the two components of the immiscible BEC. 
This task has been recently solved numerically jl5| . How- 
ever, we demonstrate that under the above mentioned 
symmetry conditions it is possible to get a reasonable 
approximation in the analytic form. First of all, let us 
represent the order parameters as 

4>i = gcosa, 4>2 = gshia. 

In the steady state, the set of Eqs.(g, g) can be than 
rewritten as 



g = —Q+ Q 

„2„/V _ 7 4 



3 ' lg^in\2a) + ga'\ (8) 



(pV) = ^ 4 sin(2a)cos(2a). 



(9) 



Here prime denotes a derivative with respect to z. The 
boundary conditions are 

a(-oo) = 0, a(+oo) = ^, q{±oo) = 1. (10) 



n=0 



,('<) 



where g( n+1 ) / g( n ) = 0(7). The second unknown func- 
tion, a, can be treated in a similar way. It is quite easy to 
find a solution in the lowest order approximation, when 
we omit the third and fourth terms in the right hand side 
of Eq.(g) and substitute g^ instead of g into Eq.(g|): 
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v^ 2 



(11) 



The "slowly varying co-ordinate" ^jz is treated here as a 
quantity of order of unity. Then we find from Eq. m) the 
first order correction to g. If we substitute g = g(°' + g^ 
into Eq.(|8j) and linearize this equation with respect to 
pW we get 
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37cos 2 a (0) sin 2 a (0) . 



Then it is important to note that gW " is of order of j 2 
rather then of 7 1 and, hence, have to be omitted when 
we find the first order correction which, finally, can be 
written as 
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Note, that only zeroth order approximation for a is nec- 
essary for determination of g^\ 

After Eq.(|l2|) has been obtained, we can find the first 
order correction to a. We approximate the latter quan- 
tity by the following formula: 



a = arctan e 



V^z+giV^z) 



(13) 



where g = 0(j). Substituting this approximation into 
Eq.(g) where g^ + (f 1 ' stands for g, after some tedious 
calculations we obtain 



37 
9= — tanh(V7^)- 



(14) 



Thus the steady state order parameter for the 1st com- 
ponent of the BEC with the O (-f 2 ) accuracy reads as 
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where 



3 7 



C = Viz + — tanh(v^z). 



(15) 
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To find 02 (z) , one simply have to apply Eq.([7|). The 
analytic approximation given by Eqs.rtlq, Hq) is justified 
by numerical calculations. The difference between the 
analytic and numerical solutions behaves as C(7 e ) in a 
wide range of 7 less than 0.1. 



To investigate the long-wavelength oscillations, we ap- 
ply the variational method analogous to that used in Ref. 
Jl6[ to study the surface oscillations of a dense BEC in 
a trap. This method implies minimization of the action 
functional associated with the Lagrangian 
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For long-wavelength perturbations, we represent the test 
functions as 



$, 



<f>j{z) + (3j(t) cos(kx)hj(z) — c/)j(z) 



ex.p[i£j(t)fj(z) cos(kx)], 



d 
' dz 



(18) 



j = 1,2. The quantity — j3j(t)hj(z) cos(kx) can be inter- 
preted as a displacement in the ^-direction of an element 
of a quantum gas that in the stationary case has the den- 
sity 4> 2 j{z). The gradient of the phase £j(t)fj(z) cos(kx) 
multiplied by two gives the velocity of the oscillating 
quantum gas. 

Due to the problem symmetry, 



h\(z) = h 2 (-z). 
The continuity equation implies that 



(19) 



hj(*)fe<t<j(*) 



dz 
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j = 1, 2. It is easy to see that it follows from Eqs.(|l9|, |20| ) 
that 



fi(z) = -M-z). 

In the long-wavelength case, when 



(21) 
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the function hj (z) varies on a spatial scale of order of k _1 . 
From the other hand, since hj(z) appears in Eq.(|l8|) only 
multiplied by dip 2 (z) j dz , one needs to know the exact be- 
haviour of hj(z) only for \z\ < 7 -1 ^ 2 , provided that hj(z) 
is finite outside this small range of z. It means that, 
because of Eq. ( |22| ) , one can choose 



h j (z) = l. 



(23) 



To find the quantity fj(z), we use the technique of 
tailoring asymptotics. Namely, for large negative z (let 
say, for z < — 7 -1 ^ 2 ), where <fi 2 (z) « 1, Eq.(E0h is reduced 
to 

^f 1 (z) = k 2 f 1 (z) 7 



and its solution approaching zero when z — > —00 is 
fi(z) — C\e kz . In the opposite case, in the range of 
z where <j)\ (z) rapidly decreases, the second term in the 
right hand side of Eq(|20|) is negligible compared to the 
first one, because of the condition Eq.(E3), and, hence, 
fi(z) = C<2 — z/2. The matching condition for these two 
asymptotics and their first derivatives allows us to deter- 
mine the constants C\ and C%. Thus, in the subsequent 
calculations we use the following expression 



fi(z) = -—e kz Q(-z) - (— 

J n ' 2k y > \2k 



+ 7: 0(^), (24) 



where 0(z) is the Heaviside's step function. This expes- 
sion and its first derivative with respect to z are continu- 
ous everywhere, including the point z = 0. The result of 
Eq.(p4h is also jus tified by numerical calculations; f%{z) 
is given by Eq.(pi|). 

Then we write the Lagrange equations for the general- 
ized coordinates (3j (t) , £j (£) . First, two of these equations 
enable us to exclude the two variables £j = $j . After this 
has been done, we arrive at the two equations containing 
f3j, [3j and describing two coupled harmonic oscillators. 
The normal modes correspond to the two cases, 



/3i=?/3 2 , 



? = ±1. 



For k <C y^ their eigenfrequencies uo q are given by 



uj* = 2 



dz {2(1 + */)(l-<;)P(z)P(-z) + 



k 2 



dz 
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}, 



dzfi(z)—(f>l(z) 
dz 



, (25) 



where P(z) — d<f>\(z)/dz. 

If the exact ground state wave function is used in 
Eq. ( p5[ ) , then lo^ is a rigorous upper bound to the oscilla- 
tion frequency for the normal modes with no nodes along 
the z-axis. Expanding the right hand side of Eq.([25|) in 
series in 7 as well as in kj -1 ' 2 and retaining the leading 
terms, we obtain the following expressions for the lowest 
oscillation frequencies for two weakly segregated BECs 
in the long wavelength regime: 



Vlk 3 , 



U-l = g\/7^- 



(26) 

(27) 



Now we can clarify the physical meaning of the two 
branches of the dispersion law Eq.(25). The q = +1 
branch corresponds to in-phase oscillations of the two 
components and represents interface bending. The wave 
number dependence, w+i ex fc 3 ' 2 for A; — s- 0, is specific 
for this type of surface waves (capillary waves) |l7J] . Also 
these interface oscillations are an analog of soft modes 
p8| in a mixture of miscible BECs with 7 being very 
small by the absolute value. We call these oscillations 
"in-phase" because for them the z-projections of the ve- 
locities of both the two components have the same sign. 



It is necessary to note that the value of A; w^ that is, in 
fact, the ratio of the surface tension constant to the sum 
of the bulk densities of the 1st and 2nd BEC components 
is equal to «/7, i.e., displays the correct dependence on 
7 following from the surface tension energy estimation 
given by Timmermans ||. 

The <; = — 1 branch corresponds to the out-of-phase os- 
cillations and, hence, to a time dependence of the depth 
of interpenetration of the two components. The disper- 
sion law, uj + 1 oc k, is analogous to that of surface waves 
on a boundary between a single component BEC and 
vacuum [|6[. 

Now it is necessary to explain why in the present prob- 
lem one cannot choose 



ui(z) = —qu-2{—z). 



(31) 
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kz 
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as it is done in Ref. [|T(|. The reason is that if one ap- 
plies Eq.(28) then for z — ► +oo the difference between 
the steady-state order parameter <f>\ (z) and the time- 
dependent order parameter < i ) i(x, z, i) is much larger (by 
the exponential factor e kz ) than <j>i{z). So, the condition 
of the relative smallness of the order parameter pertur- 

roach 
ll, the 



bation breaks down, and the validity of the ap 
described above cannot be guaranteed. In Ref. [1_^ 
condensate is assumed to fill only half of the space, and 
Eqs.(p!8|) are used for z < only, while the condensate 
density for z > is identically zero, and no such a diffi- 
culty occurs. 

In our case, the wrong choice of hj{z) and fj(z) given 
by Eq.(p8|) results in the correct value for lo 2 l 1 for k <C y^y 
but substantially, by a factor of order of 7 _1 , overesti- 
mates the eigenfrequency for the in-phase mode. The 
latter result is apparently inconsistent with the surface 
tension estimations given in Ref. ||. Moreover, use of 
Eq. (p8[) , instead of Eqs.(p3|, p4J), leads to an erroneous 
conclusion about non-monotonous behaviour of the dis- 
persion laws for both the in-phase and out-of-phase os- 
cillation when k approaches y^y. Certainly, the latter is 
simply an artefact of ineligible choice of the test function. 

Another important limiting case is the case of very 
short wavelengths, k ^> 1. Under this condition an ele- 
mentary excitation is very close to a single particle excita- 
tion (in other words, the v coefficient of the Bogoliubov's 
transformation |19| ] is small compared to the u coeffi- 
cient). Thus, the oscillation frequency is the eigenvalue 
of the Hamiltonian problem 

_<(*) + [2#(*) + (l+ 7 )#(-z)-l+ 

fc 2 ]u 1 (z) + (l+ 7 )(/> 1 (20<M-z)u 2 (z) =w«i(«), (29) 



«£(*) + [2#(-*) + (l+ 7 )#(*)-l+ 

k 2 ] u 2 (z) + (1 + 1 )cj) 1 {z)(j> 1 {-z)u 1 {z) 



um 2 (z), (30) 



It is obvious that there are two types of solutions of this 
set of coupled Scrodinger equations corresponding to the 
in-phase (q = +1) and out-of-phase (<j = —1) oscillations: 



The set of Eqs.(|29|, p0|) was solved taking into ac- 
count the symmetry condition Eq.(p^) by the Hamilto- 
nian diagonalization in the basis of harmonic functions 
L -1 / 2 sm[nm(z + L)/(2L)], where to is a positive inte- 
ger number, and L is chosen to be much larger than 
7~V 2 . rp^ e S p ec ^ rurn of the in-phase oscillations con- 
sists of the single discrete value equal to k 2 + 0.327 and 
the two branches of the continuous spectrum with the 
dispersion laws k 2 + 1 + Q 2 , Q > 0, and k 2 + 7 + q 2 , 
q > 0. One branch corresponds to the incidence to the 
intercomponent boundary of sound waves from both the 
left and right sides, and Q is the absolute value of the 
z-component of the wave vector of the incident wave. 
No atoms of one component are injected into the bulk 
of the condensate composed of atoms of another kind. 
In contrast, the second branch corresponds to injection 
of atoms of the 1st kind into the bulk of the conden- 
sate of atoms of the 2nd kind, and vice versa; q is the 
z-projection of the momentum of an injected atom. 
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FIG. 1. Numerical solutions of Eqs.(G9, KOj) the in-phase 
(solid line) and out-of-phase (dashed line) oscillations; 
7 = 0.01; z-coordinate is dimensionless, the eigenfunctions 
are plotted in arbitrary units. 

The spectrum of the out-of-phase oscillations contains 
again single discrete value equal to k 2 + 1 — 0.0447 and 
the two branches of continuous spectrum. The dispersion 
laws and physical meanings of these two branches are the 
same as for similar branches in the case of in-phase oscil- 
lations. The fact that the discrete eigenvalue correspond- 
ing to the nodeless eigenfunction shown by the dotted line 
in Fig. 1 is not the minimum eigenvalue in this case looks 
surprising, nevertheless, it is quite natural. Indeed, if we 
substitute Eq. (p^) into Eq. (E9T) we get a problem with the 
Hamiltonian containing the inversion operator. And for 
a Hamiltonian of such a kind the usual theorem, saying 
that a nodeless eigenfunction corresponds to the lowest 
eigenvalue, does not apply, as can be seen from various 
exactly solvable examples. 



The eigenfunctions u q {z) corresponding to the discrete 
eigenfrequencies are shown in Fig.l. Note, that the out- 
of-phase oscillation is much localized than the in-phase 
one. 

To summarize, we performed an analysis of dispersion 
of waves on a boundary between two weakly segregated 
dilute atomic BEC. For very small wavenumbers k we re- 
cover usual capillary and surface waves, with frequency 
squared proportional to k 3 and k, respectively. For large 
wavenumbers, each of the two types of oscillations (the 
in-phase and out-of-phase modes) has only one particu- 
lar mode well localized in z-direction near the boundary. 
Other surface modes are associated to simultaneous emis- 
sion of either sound into the bulk condensate or to free 
particles of one kind to the BEC composed of atoms of 
another kind. The latter fact must make (due to the ef- 
fect of quantum depletion) the depth of overlapping of 
the two weakly segregated ultracold bosonic clouds even 
larger than the simple mean-field theory predicts. This 
may be a very important cause of difficulties in an exper- 
imental determination of the parameter 7 fll2| ] . 
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